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I.  INTRODUCTION 

For  many  analysis  problems  in  turbomachinery,  an  assumption 
of  inviscid  flow  provides  sufficiently  accurate  results  for 
design  or  development  tasks.   At  the  same  time  the  assumption 
or  irrotational  flow  will  limit  application  to  low  speed  shock- 
less  cases.   The  Euler  equations  describe  inviscid  rotational 
flow  and  therefore  allow  the  modeling  of  transonic  or  super- 
sonic flows  with  shock  waves  and  rotational  flow  fields  behind 
strong  waves.   Solutions  of  the  Euler  equations  are  therefore 
important  for  many  applications. 

A  computer  program  was  developed  which  solves  numerically 
two-dimensional  Euler  equations  using  the  Godunov  method 
(Ref .  1)  or  Higher  Order  Extension  of  the  Godunov  method 
(Ref .  2) .   The  program  was  intended  primarily  for  turbomachinery 
applications,  but  with  modification  of  the  boundary  conditions, 
application  to  any  problem  which  is  modeled  by  the  Euler  equa- 
tions is  possible. 

The  purpose  of  this  report  is  to  give  instructions  to 
potential  users  of  the  present  computer  code  and  to  outline 
the  basic  assumptions  made  in  the  mathematical  model  and 
numerical  solution. 


II.   GOVERNING  EQUATIONS 

The  unsteady  two-dimensional  Euler  equations  can 
be  written  in  conservation  law  form  as 


8a    8b~    8c  _  Q 

at   ax   9y 


(i) 


where 
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and 
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,          c    = 
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2     2 

U   +  V 
e  =  p(e  +  ^ )  =  energy  of  the  unit  of  volume 


e  =  -r    ,  > —  _  internal  energy 

p  =  density 

U  and  V  are  the  velocity  components  corresponding  to 
the  x  and  y  coordinates 

P  =  pressure 

Y  =  specific  heat  ratio 

Equation  (1)  is  non-dimensionalized  using  a  refer- 
ence length  p  ,  a  reference  pressure x   and  a  reference 
density  p  .   These  reference  quantities  are  in  princi- 
ple arbitrary  and  were  selected  as  follows 

x   =  relative  length  of  the  geometrical  cord 

p  n   =  101350.0  -It 
m 


p   =  1.2  kg/nr 


The  reference  quantities  for  time,  t   and  velocity,  a   are 
not  independent  and  are  determined  from 


at,       p^ 

o   o  =  1  ;    __£_-_ 

x  p^  a  2 

o  o   o 


=  1  ;  (3) 


then 


Po 


t   =  x   /a 
o     o  /  o  ; 

non-dimensionalize  the  variables  in  Eq.  (1)  as  follows; 

U  =  a  U   ;  p  =  p  p   ; 
o  n  *        ^o  *n 


p  =  Pq  pn  ;  V  =  aQ  Vn 

x  =  xo  xn  '   y    =   xo  ^n 
o   n 


(4) 


When  these  quantities  are  inserted  in  Eq.  (1) ,  the  flow 
equations,  in  terms  of  non-dimensional  variables,  are 
identical  in  form  to  those  in  terms  of  dimensional 
variables.   The  subscript  "n"  may  then  be  dropped  without 
confusion  and  only  non-dimensional  quantities  need  be 
referred  to  in  the  remaining  sections  of  this  report. 
We  will  solve  Eq.  (1)  in  physical  space,  so  no  further 
transformation  of  the  equation  is  required. 


III.   BOUNDARY  CONDITIONS 

The  computational  domain  illustrated  in  Figure  1 
is  bounded  by  open  surfaces  (flow  may  cross  these  sur- 
faces) ,  solid  surfaces  and  inflow-outflow  boundaries. 
Formulation  and  implementation  of  boundary  conditions  at 
each  surface  are  discussed  in  this  section. 
3.1.   SOLID  WALL  BOUNDARY  CONDITION 

The  solid  wall  boundary  condition  appropriate  for 
inviscid  flow  computations  is  that  of  zero  mass  flux 
through  the  surface.   This  condition  is  difficult  to  im- 
plement uniquely  for  the  Euler  equations. 

In  the  program  developed,  the  user  has  a  choice  of 
two  types  of  boundary  conditions  on  the  solid  walls  as 
follows; 

a)  Solution  of  the  Riemann  problem  between  the 
point  of  the  domain  of  integration  and  its 
mirror  image  in  the  direction  normal  to  the 
wall  with  the  same  parameters  but,  with  an 
inverse  sign  of  velocity  component  normal  to 
the  wall.  (See  Fig.  2)  . 

b)  The  wall  pressure  is  derived  from  the  momen- 
tum equation 

^=^p!  (5) 

3n     R  KD) 

s 


where 


n  -  indicates  the  direction  normal  to  the 
surface 


U   =  velocity  parallel  to  the  surface 

R  =  radius  of  curvature  of  the  surface 
s 

If  we  add  the  flow  tangency  condition,  we  have  a 
uniquely  defined  boundary  condition  at  the  wall.   The 
solid  wall  boundary  conditions  are  imposed  at  the  cascade 
surface,  which  corresponds  to  the  segments  3-5  and  4-6 
in  the  computational  domain  shown  in  Figure  1. 

3.2.  PERIODIC  BOUNDARIES 

The  periodic  flow  boundary  conditions  are  imposed 
on  the  segments  1-3,  5-7,  2-4,  and  6-8  in  Fig.  1. 
By  virtue  of  the  periodicity  of  the  flow  with  respect  to 
y,  the  flow  parameters  over  segments  1-3  and  5-7  are  the 
same  as  those  over  2-4  and  6-8.   Hence,  the  segments  1-3, 
5-7,  2-4,  and  6-8  are  essentially  internal,  there  is  no 
need  to  apply  any  boundary  conditions  to  them.   In  the 
numerical  solution  additional  cells  are  attached  to  but 
outside  these  regions  as  shown  in  Figure  1.   The  flow 
parameters  in  the  extra  cells  are  set  equal  to  the  flow 
parameters  in  the  corresponding  cells  inside  the  computa- 
tional domain  which  are  displaced  in  the  Y  direction  by 
the  cascade  blade  spacing,  S. 

3.3.  INLET  CONDITIONS 

The  U  component  of  velocity,  the  flow  angle  0,  the 
gas  stagnation  enthalpy,  H,  and  the  entropy,  S,  are  kept 
constant  at  the  inlet.   This  leads  to  a  unique  definition 
of  the  flow  parameters  at  the  inlet  boundary  as  follows 


U.   =  const;   V.  /U.   =  tan9  =  const  ; 
in  in   in 


y.  P 


in 


(Y  -i)  P, 


2       2 

U     +  VT 
in      in 
+  =  H  =  const; 
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The  explicit  formula  for  V.  ,  p.  ,  and  P.   are  therefore 
r  in'   in'       in 
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yS 
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3.4.   OUTLET  CONDITIONS 

At  the  outlet  boundary  we  define  only  pressure  p   . 
and  for  all  other  flow  parameters  apply  continuation  con- 
ditions.  This  means  that  U   . ,  V   .,  and  Pout  are  set 
equal  to  the  values  of  U,  V  and  p  one  point  ahead  of  the 
outlet  boundary. 


IV.   THE  NUMERICAL  METHOD 

In  the  program  there  is  a  choice  of  using  one  of  the 
following  numerical  methods: 

a)  First  order  accurate  Godunov  Method  (Ref .  1) 

b)  Higher  order  extension  of  the  Godunov  method  (Ref. 2) 
The  higher  order  method  is  an  extension  of  the  Godunov  method 
and  was  used  as  it  is  described  in  Ref.  2.   The  piecewise- 
linear  algorithm  was  used  rather  than  piecewise-parabolic . 

Here  we  will  describe  the  Godunov  method  as  it  is  imple- 
mented by  the  program. 

Let  us  assume  that  the  grid  is  covering  the  computational 
domain  as  is  shown  in  Fig.  2(a).   We  will  also  consider  a 
fixed  grid  with  the  indexing  of  lines  as  shown  in  Fig.  2(a) . 
The  flow  parameters  related  to  the  center  of  the  cell  will 
have  a  fractional  index:  j  +  1/2,  k  +  1/2.   The  cell  boun- 
daries will  have  one  fractional  and  one  integer  index: 
j,  k  +  1/2  or  m  +  1/2,  n.   The  parameters  at  the  time  t 
will  have  subscript  index  (i.e.,  p.    ,/2   k  ,  -.  /2 )  i-n  dis- 
tinction from  the  parameters  at  the  time  t  +  At  which  will 

.  .  .  .  ,      j  +  1/2,  k  +  1/2 

have  superscript  index:   pJ     '     '  ' 

In  this  notation  Eq.  (1)  is  approximated  to  the  first 

order  of  accuracy  by  the  following  system  of  equations  in 

finite  differences: 


a  "  ai+h,k+h 


At 


3+h,k+h 


(BJ+^k    *     AYj+^,k    "   Cj+i5/k    '     **j+%,k      + 


Bj+J5,k+i5        AYj  +  l,k+Js   "   Cj  +  l,k+h        AXj  +  l,k+Js      + 


Cj+J5,k+1        AXj+*5,k+l   "  Bj+^,k+l        AYj+3s,k+l     + 


AX.    ,  .,     -   B 


AY4    .  _,,  ) 


j,k+%  j,k+?3  j  ,k+*2  j,k+?2 
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where 


'j+h,k+h   ~   h      (Xj+l,k+l   "   Xj,k)(Yj,k+l   "   Yj+l,k} 


(Xj,k+l  "  xj+l,kJ  (Yj+l,k+l  "  Yj,k} 


AYJ,k+J2  Yj,k+1  "  Yj,k 


AXJ+^,k     "      Xj+l,k  "  Xj,k 


a   = 


P 

pU 
PV 
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,       B    = 


RU 

P+RU2 
RUV 
(E+P)U 


C    = 


RU 
RUV 
P+RV2 
(E+P)V 


Here,  and  in  the  following  text  with  capital  letters 
R,  U,  V,  P,  etc.,  we  will  note  the  parameters  calculated  at 
the  edges  of  the  zones.   To  obtain  the  value  of  these  param- 
eters we  calculate  the  Riemann  problem  with  left  and  right 
states  in  the  centers  of  adjacent  cells  (Ref.  3).   For  exam- 
ple, to  calculate  the  flow  parameters  on  the  edge  j ,k+h 
we  assume  that  the  left  state  is  located  in  j-^,k+^  and  the 
right  state  in  j+h,k+h      (see  Fig.  2) .   The  velocity  vectors 
(U,V)  .  ,  ,  ,   and  (U,V)  ..,  ,  ,,  should  be  expressed  in  terms 
of  normal  and  tangent  components  of  velocity  to  the  edge 
(j,k+l)  -  (j,k).   Then  the  Riemann  problem  is  solved  for 
following  parameters: 


Left  state:   (p,  P,  Un,  V  j-^k+as 


Right  state:  (p,  P,  UR,  V^ j+h,k*h 


(8) 


where 


U   -  the  velocity  component  normal  to  the  edge 
(j,k+l)-(j,k) 

V.  -  the  velocity  component  tangent  to  the  edge 
(j,k+l)-(j,k) 

Solution  of  this  Riemann  problem  will  give: 


<R-  p<  V  Vj,k+!5 


Then  we  transform  the  component  of  the  velocity  back  to  the 
cartesian  coordinate  system  and  will  have  the  final  value  of 
the  parameters  on  the  cell  edge: 

9 


<R' p'  u<  B)j,k^ 

That  is  a  brief  description  of  the  Godunov  method. 
A  more  detailed  description  of  the  method  can  be  found  in 
Reference  4. 

The  program  includes  higher  order  extension  of  the  Godunov 
method  proposed  by  Collela  (Ref.  6)  which,  in  order  to  im- 
prove accuracy,  adds  the  following  steps  before  calculating 
the  Riemann  problem: 

1)  Compute  linear  profiles  of  the  dependent  variables 
in  each  zone  by  interpolating  slopes  at  the  centers 
of  the  zones,  subject  to  certain  monotonicity 
constraints.   This  gives  rise  to  a  global  distri- 
bution of  the  dependent  variables  which  is  piece- 
wise  linear,  linear  in  each  zone,  with  jump  dis- 
continuities at  the  edges  of  the  zones. 

2)  Correct  left  and  right  states  by  tracing  approxi- 
mate characteristics  and  solving  difference  ap- 
proximations to  the  characteristic  equations. 
Only  then  is  the  Riemann  problem  solved  for  the 
corrected  right  and  left  states  and  the  values 
obtained  by  this  solution  are  assigned  to  the 
corresponding  edge  of  the  grid  cell. 
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V.  COMPUTER  PROGRAM 

The  program  described  here  was  developed  on  the  NPS 
IBM  370-3033  computer. 

Since  the  computer  code  was  first  developed  in  Pascal, 
which  implies  structured  programming,  the  program  structure 
is  simple.   The  Data  Flow  Diagram  in  Fig.  3  illustrates  the 
general  flow  of  data  in  the  program  and  the  subroutine 
functions.   The  diagram  is  self-explanatory. 
5.1.   GENERAL  DESCRIPTION 

The  computer  program  consists  of  18  routines.   Here  we 
will  describe,  in  short,  the  general  function  of  each  routine. 
To  understand  fully  the  place  of  each  routine  and  the  program 
structure,  the  reader  should  follow  the  Data  Flow  Diagram 
in  Fig.  3  while  reading  the  written  description. 

The  subroutine  MAIN  calls  all  other  subroutines  in  the  pro- 
gram.  It  controls  the  execution  of  the  program  by  checking 
two  constants  ITRNUM  and  ITRCH.   ITRNUM  is  the  number  of 
iterations  performed.   ITRCH  controls  the  number  of  iterations 
from  the  last  call  to  the  write  routines  (WRVELC,  WRPRES, 
WRDENS  or  WRINTG) .   Three  important  constants  are  defined 
in  MAIN: 

MAX  =  number  of  grid  points  in  x  direction 

MAY  =  number  of  grid  points  in  y  direction 

LEV  =  intergration  level.   If  LEV=10  the  integration 

will  be  on  two  grid  levels  (not  recommended) .   If 
LEV=9  it  is  only  one  level  of  integration. 

Subroutine  RDDATA  reads  the  grid  definition  data  which  is 

created  on  File  3  by  the  successive  execution  of  the  programs 
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GRID  and  TRNEW.   The  program  GRID  creates  a  two  dimensional 
grid  of  the  H-type  for  the  given  cascade  definition.   The 
program  TRNEW  redefines  the  grid  to  the  form  needed  for 
application  of  the  method  and  calculates  the  cosine  and  sine 
of  angles  defining  the  direction  normal  to  the  cell  edge. 
Program  TRNEW  also  defines  the  array  REG  which  simulates  the 
computational  domain  by  assigning  to  each  physical  point  in 
the  domain  an  integer  number.   The  description  of  the  programs 
GRID  and  TRNEW  will  not  be  given  in  the  present  manual.   Here 
we  shall  assume  that  File  3,  for  a  given  cascade,  has  been  created. 
This  file  contains  the  following  data: 

REG  -  array  of  integer  numbers  which  simulates  the 
structure  of  the  computational  domain 

Al   -  array  containing  Ay  values  on  the  cell  edges  and 

the  area  of  the  cell  located  in  the  point  (j+h,k+%) 

Bl   -  array  containing  Ax  values  on  the  cell  edges 

TCOS  -  cosine  of  the  angle  between  the  normal  direction  to 
the  cell  edge  and  direction  of  the  x-axis 

TSIN  -  sire  of  the  angle  between  the  normal  direction  to 
the  cell  edge  and  direction  of  the  x-axis 

VF   -  x  coordinates  of  the  points  on  the  grid  covering 
the  computational  domain 

VY   -  y  coordinates  of  the  points  on  the  grid 

If  the  program  continues  integration  from  a  certain  step 

for  which  the  values  of  the  flow  parameters  are  known,  these 

values  should  be  stored  on  File  4  (program  creates  new  data  on 

File  4  after  a  certain  number  of  iterations) .   In  this  case 
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subroutine  RDDATA  reads  three  additional  arrays  from  File  4: 

RO  -  density 

P   -  pressure 

UB  -  velocity  component  in  x  direction 

VB  -  velocity  component  in  y  direction 
Subroutine  RDDATA  defines  the  following  constants: 

ITNO  -  number  of  iteration  in  the  solution  of  the 
Riemann  problem 

RZ  -  characteristic  density 

PZ  -  characteristic  pressure 

XZ  -  characteristic  length 

G  -  gamma 

Subroutine  DEFINE  defines  the  following  parameters : 

ITRNUM  -  number  of  iteration 

TIMTOT  -  total  time 

RADL   -  array,  radii  of  curvature  on  the  lower  side  of 
the  profile 

RADUP   -  array,  radii  of  curvature  on  the  upper  side  of 
the  profile 

If  the  program  is  executing  from  the  beginning,  the  intial  values 

of  the  array  RO,  P,  UB,  VB  are  defined  in  the  subroutine  DEFINE. 

The  initial  enthalpy  -  ENTALP,  entropy  -  ENTROP  and  tan0  -  TAU 

are  calculated  in  this  subroutine. 

Subroutine  TIME  calculates  the  time  step  for  the  next 

integration  based  on  CFL  criterion  -TJ.   This  subroutine  writes 

the  values  of  ITRNUM,  TJ,  TIMTOT,  MOMENT,  ans  MASS  where 

MOMENT  -  integral  of  momentum  in  x  direction 
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MASS  -  integral  of  mass 

Subroutine  INOUT  defines  the  inlet  and  outlet  boundary 

conditions  for  the  flow  parameters.   In  this  routine  we  cal- 

f 

culate  the  values  of  the  flow  parameter  in  the  points  (2,k) 
and  (MAX-l,k) ,  where  k  has  the  values  corresponding  only  to 
the  centers  of  the  cells. 

Subroutine  BOND  defines  boundary  conditions  for  the  com- 
putational domain,  except  inlet  and  outlet  boundary  conditions 
which  are  defined  in  the  routine  INOUT.   Two  kinds  of  boun- 
dary conditions  are  calculated  in  this  routine,  as  a  function 
of  the  cell  type  indicator  REG(j,k): 

1)  Boundary  condition  at  the  upper  and  lower  surfaces 
of  the  cascade  profile 

2)  Symmetry  boundary  conditions  for  the  regions  1-3, 
5-7,  2-4,  and  6-8.   To  calculate  the  values  of  the 
flow  parameters  at  the  cell  edges  this  subroutine 
calls  routine  ITER  to  solve  the  Riemann  problem  at 
the  edges. 

Subroutine  RIMSOL  defines  and  solves  the  Riemann  problem 
along  k  lines  between  the  centers  of  the  cells.   This  routine 
has  an  option  to  calculate  the  more  exact  values  for  the  left 
and  right  states  calling  routine  CHARY.   A  condition  state- 
ment in  this  routine  controls  call  of  the  routine  CHARY  that 
allows  the  use  of  the  same  code  for  the  Godunov  method  and  its 
higher  order  accurate  extension. 

Subroutine  RIMSLX  defines  and  solves  the  Riemann  problem 
at  the  edges  of  the  cells  in  the  y  line  directions.   In  the 
case  of  higher  order  accurate  integration,  this  routine  should 
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call  the  subroutine  CHAR.   That  call  is  controlled  by  a  con- 
dition statement. 

Subroutine  INTGR  performs  integration  of  Eq.  (.1)  in  accord- 
ance with  the  numerical  formula  (7) .  This  subroutine  calculates 
the  values  of  MOMENT,  MASS  and  ENERGY. 

Subroutine  WRVELC  calculates  and  writes  the  values  of 
the  local  Mach  number  for  the  center  of  each  cell  on  File  6. 

Subroutine  WRPRES  writes  the  values  of  the  normalized 
pressure  P(j,k)  in  the  centers  of  the  cells  on  File  6.   As 
an  option  this  subroutine  calculates  and  writes  the  pressure 
coefficient  C   instead  of  P(j,k). 

Subroutine  WRDENS  is  included  as  optional  and  it  writes 
the  values  of  the  array  RO(j,k)  (normalized  density)  on 
File  6. 

Subroutine  UUE  writes  the  values  of  the  local  Mach  num- 
ber corresponding  to  the  profile  surface  on  File  6. 

Subroutine  WRINTG  writes  complete  content  of  the  arrays? 
RO,  P,  UB  AND  VB,  on  File  4.   That  allows  in  the  next  execu- 
tion of  the  program,  beginning  integration  from  the  last  it- 
eration of  the  previous  execution  of  the  program. 

Subroutine  VISCO  is  optional  and  is  not  used  in  regular 
computation.   This  routine  adds  the  numerical  viscosity  of 
Lapidus  type  (Ref.  5)  to  the  conservation  equations . 

Subroutine  RIEMAN  calculates  the  Riemann  problem  for  the 
defined  left  and  right  states. 
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Subroutine  ITER  (called  only  from  subroutine  RIEMAN) 
calculates  pressure  on  the  interface  for  the  solution  of  the 
Riemann  problem. 

Subroutine  CHAR  corrects  values  of  the  left  and  right 
states  solving  the  equations  of  gas  dynamics  in  the  charac- 
teristics form.   In  CHAR  this  correction  is  done  in  the 
x  direction. 

Subroutine  CHARY  is  the  same  as  CHAR,  only  for  y  direction 

Subroutine  MONO,  is  called  only  form  CHAR  or  CHARY.   This 
subroutine  computes  linear  profiles  of  the  dependent  variables 
in  each  zone  by  interpolating  slopes  at  the  centers  of  the 
zones,  subject  to  certain  monotonicity  constraints. 
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VI.   INPUT 

All  external  input  to  the  program  is  written  on  File  3. 
Data  on  File  3  is  written  in  Format  (2613)  for  array  REG(j,k) 
and  in  Format  (5F15.12)  for  all  other  arrays.   If  the  user 
wants  to  continue  iteration  from  the  last  iteration  of  the 
previous  execution  of  the  program,  then  the  user  should  input 
the  values  of  the  arrays  of  the  flow  parameter.   This  data 
is  recorded  on  File  4  with  Format  (5F15.12). 

All  external  input  arrays  are  read  from  the  subroutine 
RDDATA.   The  inlet  parameters  of  the  flow  and  the  outlet  pres- 
sure are  defined  in  the  subroutine  DEFINE  from  line  1560  to 
1710. 
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VII.   OUTPUT 

All  output  of  the  program  are  recorded  on  File  6  and 
File  4.   The  following  subroutines  write  output  on  to 
File  6: 

TIME,  WRVELC,  WRPRES,  WRDENS ,  UUB . 
Only  subroutine  WRINTG  writes  on  File  4. 
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VIII.   RESULTS 

The  Fortran  version  of  the  computer  code  was  tested  for 
a  number  of  problems  for  which  the  solution  was  obtained  with 
the  Pascal  version.   The  conclusions  from  this  testing  are: 

1)  Fortran  program  is  -  50%  faster  than  Pascal  program, 
when  executed  with  the  option  Optimize  (2)  on  the 
Fortran  H  compiler  with  single  precision. 

2)  Execution  with  Fortran  G  compiler  is  30%  faster  than 
Pascal. 

3)  Execution  with  Fortran  H  compiler  with  Optimize  (1) 
option  is  30%  faster  than  Pascal. 

4)  Execution  with  Fortran  H  compiler  with  double  pre- 
cision with  Optimize  (2)  option  is  only  25%  faster 
than  Pascal. 

5)  The  differences  in  results  are  caused  by  different 
levels  of  precision  involved  in  the  calculations. 
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Fig.  1. 


Computational  grid  and  boundary  of  th< 
computational  domain. 


We  look  for  the  solution  of  the_  system  of  equations 
represented  by  equation  (1)  in  the  computational 
domain  shown  in  Fig.  1  for  t++»  with  the 
following  conditions  at  the  domain  boundaries: 

a)  Inflow  along  segment  1-2 

b)  Outflow  along  segment  7-8 

c)  Solid  wall  along  segments  3-5  and  4-6 
Cascade:  Periodicity  between  segments 

2-4  and  between  segments  5-7  and 


and 


1-3 
6-8 
Channel : 
and  2-8 


Solid  wall  along  segments  1-7 


21 


a)  Grid  Location  in  the  Coordinate  System 


K  +  1  - 


i    i+i    x 


b)  Single  cell 
Fig.  2.  Computational  Space  Notation 
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Fig.  3.   Data  Flow  Diagram  (DFD) 
of  the  program. 
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